# Section 1 工作区设置和R包加载 ----------------------------------------------------


setwd("E://OneDrive//研究//7-增温和光质光强//2-论文//1-光质与光强的关系//1-数据与分析")

library(xlsx) # Read, Write, Format Excel 2007 and Excel 97/2000/XP/2003 Files
library(plyr) # Tools for Splitting, Applying and Combining Data
library(tidyr) # Tidy Messy Data
library(agricolae) # Statistical Procedures for Agricultural Research
library(ggplot2) # Create Elegant Data Visualisations Using the Grammar of Graphics
library(gridExtra) # Miscellaneous Functions for "Grid" Graphics
library(ggpubr) # 'ggplot2' Based Publication Ready Plots


# Section 2 数据读取与参数定义 -------------------------------------------------

rm(list=ls()) #清理环境变量

FT <- c("Ht","RL","BS","AB","UB","LN","LA","LW","Chl","Chl_N","Fm.Fo","Fv.Fo","Fv.Fm") #设置响应的功能性状

Temperature <- read.xlsx("Data_20230324.xlsx", sheetName = "Temperature" , encoding = "UTF-8") #读取温度数据
Light <- read.xlsx("Data_20230324.xlsx", sheetName = "Light" , encoding = "UTF-8") #读取光照数据
BL_2 <- read.xlsx("Data_20230324.xlsx", sheetName = "Biomass_Leaf_20221031" , encoding = "UTF-8") #生物量与叶片数据
G_2 <- read.xlsx("Data_20230324.xlsx", sheetName = "Grow_20221031" , encoding = "UTF-8") #第2次生长数据
G_1 <- read.xlsx("Data_20230324.xlsx", sheetName = "Grow_20221013" , encoding = "UTF-8") #第1次生长数据

Data_last <- merge(G_2, BL_2, by = c("Num", "Pot", "C", "T", "L", "rep", "Num_plant", "S")) #收苗时得到的全部数据
Grow_twice <- cbind(data.frame(Date = c(rep(1013, 300), rep(1031, 300))), rbind(G_1, G_2)) #两次生长数据

# Section 3 自定义函数 ---------------------------------------------------------

se.yb <- function(c) {
  c <- c[!is.na(c)]
  se <- sd(c)/sqrt(length(c))
  return(se)
} #标准误计算
  
temperature.yb <- function(data, color) {
    library(ggplot2) # Create Elegant Data Visualisations Using the Grammar of Graphics
    library(gridExtra) # Miscellaneous Functions for "Grid" Graphics
    library(ggpubr) # 'ggplot2' Based Publication Ready Plots
    data$d <- data$T - data$CK
    Fig_Tsum_data <- data.frame(A = c("Ambient temperature(°C)", "Warming temperature(°C)", 
                                      "Temperature difference(°C)", "T-value", "Significance"), 
                                B = c(paste(round(mean(data$CK),2), "±", 
                                            round(sd(data$CK)/sqrt(length(data$CK)),2)),
                                      paste(round(mean(data$T),2), "±", 
                                            round(sd(data$T)/sqrt(length(data$T)),2)),
                                      paste(round(mean(data$d),2), "±", 
                                            round(sd(data$d)/sqrt(length(data$d)),2)),
                                      round(t.test(data$CK, data$T, pair=TRUE)[["statistic"]],3),
                                      round(t.test(data$CK, data$T, pair=TRUE)[["p.value"]],3)))
    if (Fig_Tsum_data[5,2] < 0.001) {
      Fig_Tsum_data[5,2] <- "<0.001"
    }
    g_data <- data.frame(date = rep(data$date,2),
                         temperature = c(data$CK, data$T),
                         group = c(rep("Ambient", length(data$date)), 
                                   rep("Warming", length(data$date))))
    Fig_T <- ggplot(g_data)+
      geom_line(aes(x = date, y = temperature, color = group), linewidth = 1)+
      scale_x_date(date_breaks = "1 week",
                   date_labels = "%m/%d")+
      scale_y_continuous(limits = c(min(g_data$temperature), 
                                    1.25*max(g_data$temperature) - 0.25*min(g_data$temperature)))+
      labs(x = "", y = "Temperature(°C)")+
      scale_color_manual(values = color)+
      annotation_custom(grob = tableGrob(t(Fig_Tsum_data), rows = NULL, cols = NULL),
                        ymin = 1.1*max(g_data$temperature) - 0.1*min(g_data$temperature),
                        ymax = 1.2*max(g_data$temperature) - 0.2*min(g_data$temperature))+
      theme_bw()+
      theme(legend.position="top",
            legend.title = element_blank(),
            legend.text = element_text(size = 12, face = "bold"),
            axis.text.x = element_text(size = 12, face = "bold"),
            axis.title.x = element_blank(),
            axis.title.y = element_text(size = 14),
            axis.text.y = element_text(size = 12))
    return(Fig_T)
  } #温度数据处理

transpose.yb <- function(data, col_DV, col_IV) {
  for (i in c(1:length(col_DV))) {
    FT = rep(col_DV[i], length(data[,1]))
    df <- cbind(FT, data[,c(col_IV, col_DV[i])])
    names(df) <- c("FT", col_IV, "DV")
    if (i == 1) {
      out <- df
    } else {
      out <- rbind(out, df)
    }
  }
  return(out)
} #转为单列数据，原数据头做一新列
  
sensitivity.yb <- function(data, col_Treat, value_Treat0, col_IV, col_DV, col_group) {
  library(tidyr) # Tidy Messy Data # Tidy Messy Data
  library(plyr) # Tools for Splitting, Applying and Combining Data
  df <- data
  
  if (!is.null(col_group)) {
    df <- unite(df, "group", all_of(col_group), sep = "_", remove = F)
    df$group <- as.factor(df$group)
    group <- unique(df$group)
  } else {
    df$group <- as.factor(1)
    group <- unique(df$group)
  } #合并确定分组
  
  #计算dx
  for (i in c(1:length(col_Treat))) {
    if (i == 1) {
      dx <- df[,col_IV[i]] / unique(df[df[,col_Treat[i]] == value_Treat0[i],col_IV[i]])
    } else {
      dx <- dx * df[,col_IV[i]] / unique(df[df[,col_Treat[i]] == value_Treat0[i],col_IV[i]])
    }
  }
  df$dx <- dx-1
  
  
  #计算dy
  data_list <- list()
  out_list <- list()
  Data <- data.frame()
  for (i in c(1:length(group))) {
    data_list[[i]] <- df[df$group == group[i],]
    out_list[[i]] <- df[df$group == group[i], c(col_Treat, "dx", col_group)]
    for (j in c(1:length(col_Treat))) {
      if(j == 1) {
        CK_row <- which(data_list[[i]][,col_Treat[j]] == value_Treat0[j])
      } else {
        CK_row <- intersect(CK_row, which(data_list[[i]][,col_Treat[j]] == value_Treat0[j]))
      }
    }
    for (j in c(1:length(col_DV))) {
      y0 <- mean(data_list[[i]][CK_row,col_DV[j]])
      dy <- (data_list[[i]][,col_DV[j]] - y0)/y0 *100
      out_list[[i]] <- cbind(out_list[[i]], dy)
      names(out_list[[i]])[names(out_list[[i]]) == "dy"] <- col_DV[j]
    }
    Data <- rbind(Data, out_list[[i]])
  }
  
  Data <- Data[Data$dx != 0,] #清除不必要的i=0的数据
  
  # 转置数据
  for (i in c(1:length(col_DV))) {
    FT = rep(col_DV[i], length(Data[,1]))
    Data2 <- cbind(FT, Data[,c(col_group, "dx", col_DV[i])])
    names(Data2) <- c("FT", col_group,"dx", "dy")
    if (i == 1) {
      SC <- Data2
    } else {
      SC <- rbind(SC, Data2)
    }
  }
  
  # 计算敏感度
  SC$SC <- SC$dy / SC$dx
  
  col_group <- c(col_group, "FT")
  n <- length(col_group)
  out <- list()
  out[[n+1]] <- SC
  while (n > 0) {
    out[[n]] <- ddply(out[[n+1]], eval(parse(text = paste(".(", paste(col_group[1:n], collapse = ", "), ")"))),
                      summarize, 
                      sc = mean(SC, na.rm = TRUE), 
                      N  = length(SC),
                      sd = sd(SC, na.rm = TRUE),
                      se = se.yb(SC))
    names(out[[n]])[names(out[[n]]) == "sc"] <- "SC"
    n <- n-1
  }
  out[[length(out)]] <- out[[length(out)]][, c(col_group, "SC")]
  return(out)
} #敏感度计算




# Section 4 数据处理 ----------------------------------------------------------

## 环境处理数据
### 温度
Temperature_mean <- data.frame(
  T = c("CK", "T"),
  Temperature = c(mean(Temperature$CK), mean(Temperature$T))
)

### 光照条件
#### 光强
Light_intensity <- ddply(Light, .(L), summarize,
                         Intensity = mean(E.lx.), SD = sd(E.lx.), SE = se.yb(E.lx.), N = length(E.lx.))
HSD <- HSD.test(aov(E.lx. ~ L, data = Light),"L")$group
HSD$L <- rownames(HSD)
Light_intensity <- merge(Light_intensity, HSD, by = "L")[,c("L","Intensity","SD","SE","N","groups")]

#### 光质
Light_nm <- transpose.yb(Light, paste("X", c(380:780),"nm", sep = ""), c("T","L","rep"))
Light_quality <- merge(
  ddply(Light_nm[Light_nm$FT %in% paste("X", c(650:670),"nm", sep = ""),], .(L), summarize, R  = mean(DV)),
  ddply(Light_nm[Light_nm$FT %in% paste("X", c(720:740),"nm", sep = ""),], .(L), summarize, FR = mean(DV)),
  by = "L"
)
Light_quality$R.FR = Light_quality$R / Light_quality$FR

#### 光照条件结果
light <- merge(Light_intensity[,c("L", "Intensity")], Light_quality[,c("L", "R.FR")]) 

#### 处理强度数据与性状响应值合并
Data_last <- merge(merge(Data_last, Temperature_mean, by = "T"), light, by = "L")

#### 光强组敏感度计算
SC_IT <- sensitivity.yb(Data_last[Data_last$L %in% c("N", "W" ,"S"),], 
                               "T", "CK", "Temperature", FT,  c("S", "C")) #对温度的敏感度
SC_IL <- sensitivity.yb(Data_last[Data_last$L %in% c("N", "W" ,"S"),], 
                               "L", "N", "Intensity", FT,  c("S", "C")) #对光强的敏感度
SC_ITL <- sensitivity.yb(Data_last[Data_last$L %in% c("N", "W" ,"S"),], 
                                c("T", "L"), c("CK", "N"), c("Temperature", "Intensity"), FT, c("S", "C")) #对温度和光强交互的敏感度

#### 光质组敏感度计算
SC_QT <- sensitivity.yb(Data_last[Data_last$L %in% c("N", "R" ,"RF"),], 
                               "T", "CK", "Temperature", FT,  c("S", "C")) #对温度的敏感度
SC_QL <- sensitivity.yb(Data_last[Data_last$L %in% c("N", "R" ,"RF"),], 
                               "L", "N", "Intensity", FT,  c("S", "C")) #对光质的敏感度
SC_QTL <- sensitivity.yb(Data_last[Data_last$L %in% c("N", "R" ,"RF"),], 
                                c("T", "L"), c("CK", "N"), c("Temperature", "Intensity"), FT, c("S", "C")) #对温度和光质交互的敏感度

SC_T <- merge(SC_IT[[3]], SC_QT[[3]], by = c("S", "C", "FT"))
names(SC_T) <- c("S", "C", "FT", "SC_I", "N_I", "sd_I", "se_I", "SC_Q", "N_Q", "sd_Q", "se_Q")
SC_L <- merge(SC_IL[[3]], SC_QL[[3]], by = c("S", "C", "FT"))
names(SC_L) <- c("S", "C", "FT", "SC_I", "N_I", "sd_I", "se_I", "SC_Q", "N_Q", "sd_Q", "se_Q")
SC_TL <- merge(SC_ITL[[3]], SC_QTL[[3]], by = c("S", "C", "FT"))
names(SC_TL) <- c("S", "C", "FT", "SC_I", "N_I", "sd_I", "se_I", "SC_Q", "N_Q", "sd_Q", "se_Q")

#### 相关分析

cor_T <- ddply(SC_T, .(S, C), summarise, 
               N = length(SC_I),
               R = cor.test(SC_Q, SC_I)[["estimate"]][["cor"]],
               df = cor.test(SC_Q, SC_I)[["parameter"]][["df"]],
               P = cor.test(SC_Q, SC_I)[["p.value"]])


cor_L <- ddply(SC_L, .(S, C), summarise,  
               N = length(SC_I),
               R = cor.test(SC_Q, SC_I)[["estimate"]][["cor"]],
               df = cor.test(SC_Q, SC_I)[["parameter"]][["df"]],
               P = cor.test(SC_Q, SC_I)[["p.value"]])

cor_TL <- ddply(SC_TL, .(S, C), summarise,  
                N = length(SC_I),
                R = cor.test(SC_Q, SC_I)[["estimate"]][["cor"]],
                df = cor.test(SC_Q, SC_I)[["parameter"]][["df"]],
                P = cor.test(SC_Q, SC_I)[["p.value"]])


# Section 5 绘图 ------------------------------------------------------------
SC_T <- unite(SC_T, "group", c("S", "C"), sep = "_", remove = F)
ggplot(SC_T)+
  geom_point(aes(x = SC_Q, y = SC_I, color = group))+
  geom_smooth(method = lm, aes(x = SC_Q, y = SC_I, color = group), se = F)
SC_L <- unite(SC_L, "group", c("S", "C"), sep = "_", remove = F)
ggplot(SC_L)+
  geom_point(aes(x = SC_Q, y = SC_I, color = group))+
  geom_smooth(method = lm, aes(x = SC_Q, y = SC_I, color = group), se = F)
SC_TL <- unite(SC_TL, "group", c("S", "C"), sep = "_", remove = F)
ggplot(SC_TL)+
  geom_point(aes(x = SC_Q, y = SC_I, color = group))+
  geom_smooth(method = lm, aes(x = SC_Q, y = SC_I, color = group), se = F)

# temperature.yb(Temperature, c("#1f78b4", "#33a02c"))

SC_IT

SC_T <- merge(SC_IT[[4]], SC_QT[[4]], by = c("S", "C", "FT"))
names(SC_T) <- c("S", "C", "FT", "SC_I", "SC_Q")

